% To analyze the time series related to the
% Bhuj EQ of 26-1-2001

file = 'c:\manoj\projects\eq\eq.rts';
fid = fopen(file,'rb');
[TS,er] = rts_in(fid);
ProcDef = SetDefault(TS.head.band);
ProcDef.file1=file;
%  ProcDef.block = 1024; %!!!!!! REMOVE THIS LINE 24.7.05 TOULOUSE
ProcDef.len = floor(length(TS.matrix)/ProcDef.block);
ProcDef.SelStacks = ones([1 ProcDef.len]);

for i = 1:TS.head.nchannel,
factor(i) = NrFct(TS.head.band,ProcDef.sensor(i),TS.probe.dist(i),TS.probe.gain(i));
end;
factor(1:2) = factor(1:2)*0.040; %these factors were found necessary while 
factor(3:5) = factor(3:5)*5; % comparing with ProcMT 22.5.3

waven  = 'cgau4';
nchan = 1;
%ax = [6e4:8e4];
ax = [1:length(TS.matrix)];
subplot(211);
plot(TS.matrix(ax,nchan)*factor(nchan));
axis([-inf inf -inf inf]);

%Sc = scal2frq([8 6 4 2 1 0.5 0.25 0.1],waven,1/32); %to get desired scales for a set of frequencies
Sc = scal2frq(logspace(-1,1.2,10),waven,1/32); %to get desired scales for a set of frequencies


%c = cwt(matrix(:,5),[2 4 16 32 64 128],'cgau4'); %Complex wavelet transform
c = cwt(TS.matrix(ax,nchan)*factor(nchan),Sc,waven);
%F = scal2frq([2 4 16 32 64 128 200 256],'cgau4',1/32);
F = scal2frq(Sc,waven,1/32); %should be same as the frequency set
%F = scal2frq([2 4 8 16 32 64 128],'cmor1-1',1/32);
subplot(212);
imagesc(ax,F,(log10(abs(c)))); %problem - uncalibrated ts
% Note 8 Nov 2005
%[x,y] = sinu;
%c = cwt(y(1:1000),[10,15,20.3125,25,30,35,40.6250,45,50,60],'morl');
%F = scal2frq([10,15,20.3125,25,30,35,40.6250,45,50,60],'morl',1/50);
%contour(repmat([1:1000],[length(F),1]),repmat(F',[1,1000]),c);
% with the above couple of codes, I could see that the scal2frq eturns the
% correct fft frequencies this case around 1 and 2 Hz - as given in sinu.m. 
% Only problem is the smeared power over target
% frequencies

%plotting of time series
% ax = [1:length(TS.matrix)];
% axt = ax/32;
% nchan = 1;
% subplot(511);
% plot(axt,TS.matrix(ax,nchan)*factor(nchan));
% axis([-inf inf -2 2]);
% nchan = 2;
% ylabel('mV/km');
% subplot(512);
% plot(axt,TS.matrix(ax,nchan)*factor(nchan));
% axis([-inf inf -2 2]);
% ylabel('mV/km');
% nchan = 3;
% subplot(513);
% plot(axt,TS.matrix(ax,nchan)*factor(nchan));
% axis([-inf inf -2 2]);
% ylabel('nT/s');
% nchan = 4;
% subplot(514);
% plot(axt,TS.matrix(ax,nchan)*factor(nchan));
% axis([-inf inf -2 2]);
% ylabel('nT/s');
% nchan = 5;
% subplot(515);
% plot(axt,TS.matrix(ax,nchan)*factor(nchan));
% axis([-inf inf -2 2]);
% ylabel('nT/s');
% xlabel('Time (s)');

% Plot of wavelet sepctra - 2005 version


waven  = 'cgau4';
nchan = 2;
%ax = [6e4:8e4];
ax = [1:length(TS.matrix)];
subplot(211);
plot(TS.matrix(ax,nchan)*factor(nchan));
axis([-inf inf -inf inf]);

%Sc = scal2frq([8 6 4 2 1 0.5 0.25 0.1],waven,1/32); %to get desired scales for a set of frequencies
Sc = scal2frq(logspace(-1,1.2,10),waven,1/32); %to get desired scales for a set of frequencies


%c = cwt(matrix(:,5),[2 4 16 32 64 128],'cgau4'); %Complex wavelet transform
c = cwt(TS.matrix(ax,nchan)*factor(nchan),Sc,waven);
%F = scal2frq([2 4 16 32 64 128 200 256],'cgau4',1/32);
F = scal2frq(Sc,waven,1/32); %should be same as the frequency set
%F = scal2frq([2 4 8 16 32 64 128],'cmor1-1',1/32);
subplot(212);
imagesc(ax,F,(log10(abs(c)))); %problem - uncalibrated ts
axt = ax/32;
ax = [1:length(TS.matrix)];
xdata = (log10(abs(TS.matrix(ax,nchan)*factor(nchan))));
axk = [1:1000:length(xdata)]/32;
icount=1;
for i = 1:1000:length(xdata)-1000,
k(icount) = median(xdata(i:i+999));icount=icount+1;
end;
subplot(2,1,1);
imagesc(ax/32,F,(log10(abs(c))));
 ylabel('Frequency Hz');
 xlabel('Time S');
 colorbar('horiz');
subplot(514);
plot(axk(1:end-1),k,'b');
ylabel('log(|Hz|) - nT/s');
axis([-inf 5120 -inf inf]);
%put 0.10   0.54    2.93   15.85  these
%vaues to labels for y axis (edit the image)

%plot spectra of all the channels
for i = 1:5,
waven  = 'cgau4';
nchan = i;
%ax = [6e4:8e4];
ax = [1:length(TS.matrix)];
%Sc = scal2frq([8 6 4 2 1 0.5 0.25 0.1],waven,1/32); %to get desired scales for a set of frequencies
Sc = scal2frq(logspace(-1,1.2,10),waven,1/32); %to get desired scales for a set of frequencies

%c = cwt(matrix(:,5),[2 4 16 32 64 128],'cgau4'); %Complex wavelet transform
c = cwt(TS.matrix(ax,nchan)*factor(nchan),Sc,waven);
%F = scal2frq([2 4 16 32 64 128 200 256],'cgau4',1/32);
F = scal2frq(Sc,waven,1/32); %should be same as the frequency set
%F = scal2frq([2 4 8 16 32 64 128],'cmor1-1',1/32);
axt = ax/32;
ax = [1:length(TS.matrix)];
xdata = (log10(abs(TS.matrix(ax,nchan)*factor(nchan))));
axk = [1:1000:length(xdata)]/32;
subplot(5,1,i);
imagesc(ax/32,F,(log10(abs(c(1:9,:)))));
% ylabel('Frequency Hz');
% xlabel('Time S');
caxis([-6 0.1]);
colorbar;
end;
%put 0.10   0.54    2.93   15.85  these
%vaues to labels for y axis (edit the image)